use     values-secularism.dta , replace 

* This script needs estout to cr

* create secularism variable from individual items 
* use all nonmissing items
alpha sec2-sec5
egen secularism =rowmean(sec2-sec5)

* correlation between aspects of religiosity and tradition/self-direction 
gen none=5.denomination
lab var none "Religion: none"

* make a nice table
estpost correlate tradition religin prayer church none
eststo trad 

estpost correlate selfdirection religin prayer church none
eststo self

esttab trad self using corrtable.tex ,b(2) noobs label nonumb mtitles(Tradition Self-Direction) nonotes addnotes("N=1938--1973 ") replace nogaps not alignment(S)

reg tradition i.denomination
reg selfdirection i.denomination

corr tradition prayer
corr selfdirection prayer

corr tradition churchattenda
corr selfdirection churchattenda 

corr tradition religintensity
corr selfdirection religintensity

* We can use bivariate regression and model stacking to make a nice table
* Still need to deal with the number of cases 



* Ok. run fullest regression model 

reg secularism  i.male i.east i.agecat i.educ ib5.denomination religintensity prayer churchattenda tradition selfdirect

* Get vif. Calculate max and mean and add to model 
estat vif
local maxvif = max(r(vif_1), r(vif_2), r(vif_3), r(vif_4), r(vif_5), r(vif_6), r(vif_7), r(vif_8), r(vif_9), r(vif_10), r(vif_11), r(vif_12), r(vif_13), r(vif_14))
local meanvif = (r(vif_1) + r(vif_2) + r(vif_3) + r(vif_4) + r(vif_5) + r(vif_6) + r(vif_7) + r(vif_8) + r(vif_9) + r(vif_10) + r(vif_11) + r(vif_12) + r(vif_13) + r(vif_14))/14
estadd scalar max_vif = `maxvif'
estadd scalar mean_vif = `meanvif'

* Store model and sample 
est store full
gen included = e(sample)


* run/store model with socio-demographics and religion but w/o values
reg secularism  i.male i.east i.agecat i.educ ib5.denomination religintensity prayer churchattenda if included 
estat vif
local maxvif = max(r(vif_1), r(vif_2), r(vif_3), r(vif_4), r(vif_5), r(vif_6), r(vif_7), r(vif_8), r(vif_9), r(vif_10), r(vif_11), r(vif_12))
local meanvif = (r(vif_1) + r(vif_2) + r(vif_3) + r(vif_4) + r(vif_5) + r(vif_6) + r(vif_7) + r(vif_8) + r(vif_9) + r(vif_10) + r(vif_11) + r(vif_12))/12
estadd scalar max_vif = `maxvif'
estadd scalar mean_vif = `meanvif'
est store relig

* run/store model with socio-demographics only 
reg secularism i.male i.east i.agecat i.educ if included 
estat vif
local maxvif = max(r(vif_1), r(vif_2), r(vif_3), r(vif_4), r(vif_5), r(vif_6))
local meanvif = (r(vif_1) + r(vif_2) + r(vif_3) + r(vif_4) + r(vif_5) + r(vif_6) )/6
estadd scalar max_vif = `maxvif'
estadd scalar mean_vif = `meanvif'
est store sociodem

* Make labels a bit nicer for esttab
lab var tradition Tradition
lab var selfdirection "Self-direction"
lab def east 0 "west" 1 " east",replace

* Output 
esttab sociodem relig full using linmodels.tex, nonumbers mtitles("M1" "M2" "M3") nobaselevels refcat(1.denomination "Religion:" 2.education "Education:" 2.agecat "Age group:" 1.east "Region:" 1.male "Gender:",nolabel) label stats(N r2 bic max_vif mean_vif,fmt(0 2 1 2 2) labels("N" "R2" "BIC" "VIF (max)" "VIF (mean)")) se substitute(R2 \$R^2\$) alignment(S) replace nogaps
